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1.  Introduction 


Traditional  metal  plasticity  models,  formulated  in  terms  of  strain  rates  and  stresses  and 
incorporated  in  large-scale  numerical  analyses,  provide  useful  solutions  for  a  wide  range  of 
problems.  Details  of  the  material  microstructure  interactions  that  govern  the  deformation 
response  are  assumed  to  occur  at  length  scales  not  resolvable  by  the  simulations  and  are  captured 
implicitly  in  the  constitutive  relations.  For  example,  dependence  of  the  yield  strength  on  grain 
size  through  the  Hall-Petch  effect  can  be  incorporated  by  including  grain  size  in  the  constitutive 
model  without  tracking  the  details  of  dislocation  interactions  with  grain  boundaries. 

In  simulations  with  resolution  at  the  grain  scale,  as  in  multiscale  modeling,  the  length  scales 
dictating  some  hardening  mechanisms  are  comparable  to  the  numerical  discretization.  The 
torsion  experiments  of  Fleck  and  Hutchinson  (1997)  clearly  demonstrate  increased  strength  with 
decreasing  size  for  wires  10’s  of  microns  in  diameter.  The  strengthening  is  attributed  to  gradients 
in  the  crystal  lattice  orientation,  which  create  boundaries  where  dislocations  accumulate.  These 
structures  both  store  energy  and  provide  resistance  to  further  dislocation  motion  (e.g.,  Baskaran 
et  ah,  2010).  Applying  traditional  crystal  plasticity  models  (e.g.,  Asaro,  1983;  Peirce  et  ah,  1983) 
to  investigate  the  strengthening  in  a  multiscale  framework  will  not  be  successful  because  the 
models  are  formulated  in  terms  of  traditional  continuum  variables  of  strain  rate  and  stress,  and 
there  is  no  underlying  microstructure  length  scale  that  would  produce  a  size  effect. 

Numerous  studies  over  the  past  decade  have  investigated  ways  to  incorporate  a  length  scale  into 
the  continuum  crystal  plasticity  model.  Most  focus  on  microstructure  gradients,  as  it  is 
recognized  that  both  the  Hall-Petch  effect  and  the  results  of  Fleck  and  Hutchinson  are  tied  to 
gradients.  To  cite  just  a  few  of  the  many  examples,  models  have  examined  continuously 
distributed  dislocations  (Acharya,  2001),  dislocation  density  gradients  (Arsenlis  et  ah,  2004),  and 
gradient  related  state  variables  (Gurtin  et  ah,  2007;  Gerken  and  Dawson,  2008;  Mayeur  et  ah, 
2011).  These  formulations  all  introduce  additional  variables  into  the  solution,  which  creates  two 
difficulties.  The  first  is  fitting  these  into  a  solution  scheme  with  appropriate  evolution  equations, 
and  the  second  is  determining  an  appropriate  prescription  for  boundary  conditions  on  the  new 
variables.  The  former  is  typically  the  focus  of  the  research.  The  latter  often  remains  an  open 
issue.  For  example,  what  is  the  value  of  a  strain  gradient  at  a  free  surface  versus  an  interface,  or 
how  should  dislocation  density  be  prescribed  as  a  boundary  condition? 

When  using  the  traditional  crystal  plasticity  model  in  a  finite  element  code,  the  deformation  of 
neighboring  elements  is  only  connected  through  shared  nodal  displacements  and  force 
equilibrium  at  the  nodes.  Even  though  dislocations  associated  with  slip  travel  through  the 
material,  the  transmission  of  dislocations  from  one  element  to  another  is  not  represented,  and  the 
resulting  coordination  in  slip  deformation  is  not  captured.  The  effective  result  is  that  all  finite 
element  boundaries  are  infinite  sources  and  sinks  for  dislocations.  Advanced  crystal  plasticity 
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formulations,  such  as  those  based  on  lattice  orientation  gradients,  often  include  the  continuity  of 
dislocation  flux  across  finite  element  boundaries  as  a  byproduct.  However,  the  importance  of  this 
physical  constraint  is  not  explored  independently. 

The  goal  of  the  current  work  is  to  determine  if  a  simpler  approach  can  be  used  to  address  the 
length  scale  and  other  deficiencies  of  the  classical  crystal  plasticity  model.  It  is  a  search  in 
pursuit  of  Occam’s  razor,  the  law  of  parsimony.  The  starting  point  is  enforcing  dislocation  flux 
continuity  across  finite  element  boundaries.  This  uses  existing  solution  variables  and  the 
boundary  conditions  are  conceptually  straightforward.  Dislocation  flux  is  unconstrained  at  free 
surfaces  and  zero  at  rigid  boundaries  and  intermediate  at  grain  boundaries  that  are  sources  and 
sinks  for  dislocations.  The  flux  gradients  can  be  used  to  infer  the  evolution  of  dislocation 
gradients  that  lead  to  lattice  orientation  gradients. 

This  report  presents  the  formulation  and  results  from  the  first  year’s  effort  on  the  two-year 
Director’s  Research  Initiative  (DRI)  project.  The  finite  element  code,  the  crystal  plasticity 
model,  and  the  dislocation  flux  continuity  relations  are  presented  in  section  2.  Section  3  contains 
results  from  verification  calculations  and  initial  simulations  exploring  the  effects  of  flux 
continuity  and  hardening  associated  with  slip  gradients.  The  report  concludes  with  a  discussion 
of  the  results  and  several  issues  that  were  identified. 


2.  Finite  Element  Code,  Crystal  Plasticity  Model,  and  Slip  Continuity 


The  proposed  model  is  non-local,  where  the  material  response  in  a  finite  element  depends 
directly  on  the  structure  evolution  in  neighboring  elements.  Only  a  few  existing  finite  element 
frameworks  would  accommodate  such  a  model  without  substantial  modifications.  Another  goal, 
that  further  restricts  the  use  of  existing  software,  is  an  explicit  dynamic  implementation  that  can 
be  used  in  high  rate  deformation  applications.  Given  the  requirements,  it  was  determined  that  it 
would  be  more  efficient  and  more  effective  to  create  a  simple,  explicit  finite  element  research 
code  rather  than  to  modify  an  existing  code  and  its  data  structures.  The  finite  element  model  is 
described  in  section  2.1. 

Modules  for  the  underlying  crystal  plasticity  model  are  taken  from  an  existing  implementation  in 
ALE3D  (Becker,  2004).  A  brief  description  is  given  below.  An  idealized,  two-dimensional  (2-D) 
crystal  geometry  with  three  slip  systems,  oriented  60°  apart,  was  implemented  so  that  the 
formulation  and  algorithm  development  could  be  carried  out  more  efficiently.  Time  integration 
of  the  constitutive  model  within  the  finite  element  code  departs  significantly  from  standard  finite 
element  time  integration  schemes  and  is  described  in  section  2.4. 
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2.1  Explicit  Finite  Element  Framework 

A  2-D,  explicit  finite  element  code  was  created  as  a  framework  to  develop  and  evaluate  the  slip 
continuity  model.  It  uses  either  constant  strain  triangle  elements  in  a  crossed-triangle 
configuration  (Nagtegaal,  Parks,  and  Rice,  1974)  or  constant  strain  quadrilaterals  with  hourglass 
control  (Flanagan  and  Belytschko,  1981).  Time  integration  is  with  the  second-order  accurate 
“leap-frog”  method,  where  positions,  forces,  and  accelerations  are  evaluated  at  full  time  step 
intervals  and  velocities  are  evaluated  at  the  midpoints  of  the  time  steps.  In  the  absence  of  body 
forces,  the  momentum  equation  determines  accelerations  by 


xt  =  (y  •  o)ipt  •  (i) 

Updates  to  the  velocity  and  position  are  evaluated,  respectively,  as 

Xt+At/2  =  *t— At/2  +  XtAt  (2) 

and 


■^■t+At  Xf  . 


Written  in  the  discretized  finite  element  context,  equation  1  becomes 


r  M 


•  bjm  V 


M 


/mJ , 


(3) 

(4) 


where  B,M is  the  derivative  of  the  finite  element  shape  function  for  node  J  of  element  M,  m]  is 
the  node  mass,  and  VM  is  the  element  volume.  The  summation  is  over  all  of  the  elements 
surrounding  the  node.  Equation  4  implies  a  lumped  mass  matrix.  The  nodal  mass  contribution 
from  each  element  is  determined  from  integration  of  the  element  volume  rather  than  simply 
dividing  the  element  mass  by  four.  This  properly  accounts  for  the  element  shape  and  provides  a 
more  consistent  mass  for  the  Flanagan-Belytschko  element. 


The  velocity  gradient  for  element  M  is  calculated  using  the  mid-step  configuration  when 
determining  BJM  in  order  to  achieve  second-order  accurate  strain  integration, 


Lm  = 


(5) 


2.2  Crystal  Kinematics 

The  crystal  plasticity  model  used  in  the  current  study  follows  from  the  widely  used  kinematic 
framework  described  by  Asaro  (1983)  and  first  implemented  in  a  finite  element  framework  by 
Peirce,  Asaro,  and  Needleman  (1983).  The  deformation  gradient,  F,  is  notionally  decomposed 
into  elastic  and  plastic  parts,  Fe  and  Fp,  respectively 

F  =  Fe  ■  Fp  (6) 
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The  elastic  part  accounts  for  distortion  and  rotation  of  the  crystal  lattice,  and  the  plastic  part 
represents  slip  on  predefined  crystal  planes  and  directions  that  moves  material  but  does  not  alter 
the  underlying  crystal  lattice.  This  is  shown  schematically  in  figure  1 .  The  plastic  part  can  be 
envisioned  as  the  accumulated  effect  of  slip  on  each  slip  system,  ya,  with  the  reference 
crystallographic  plane  normals  given  by  wIq  and  reference  slip  directions  denoted  by  Sq. 

Nsys 

Fp  =  ^  ya  sa  0  ma  +  /  (7) 

a- 1 

The  elastic  part  of  the  deformation  gradient  orients  the  lattice  in  the  current  configuration  and 
distorts  it  consistent  with  the  applied  stresses. 


Figure  1.  Representation  of  the  elastic-plastic  decomposition 
of  the  deformation  gradient. 

The  intermediate  configuration  illustrated  in  figure  1  does  not  really  exist,  but  it  serves  as  a 
convenient  reference  frame  for  constructing  the  slip  model.  It  is  important  to  note  that,  while 
deformation  gradient  describes  a  compatible  deformation  field  with  no  holes  and  a  smoothly 
varying  displacement  field,  the  notional  deformation  fields  associated  with  the  elastic  and  plastic 
parts  do  not.  A  crystal  subjected  to  a  non- affine  (inhomogeneous)  deformation,  and  conceptually 
unloaded  elastically  to  the  intermediate  configuration  shown  in  figure  1 ,  may  have  cracks,  holes, 
and/or  overlaps.  In  addition,  the  total  slip  form  represented  in  equation  7  and  figure  1  is  only 
applicable  for  small  deformations  or  if  the  relative  slip  rates  on  the  various  slip  systems  are 
constant  over  the  deformation  history.  The  rate  form  of  the  model  is  more  general. 

The  velocity  gradient,  L,  resulting  from  this  kinematic  description,  creates  an  additive 
decomposition  into  elastic  and  plastic  parts: 

L  =  F  F~1  =  Fe  Fe_1  +  Fe  •  Fp  •  Fp_1  •  Fe_1  dAf  Le  +  Lp  (8) 
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Applying  the  elastic  part  of  the  velocity  gradient  indicated  by  equation  8  to  the  rate  form  of 
equation  7  gives  the  plastic  part  of  the  velocity  gradient  as 

Nsys  NSyS 

Lv  =  ^  ya(Fe  •  $o  ®  •  Fe_1)  =  ma ,  (9) 

a-1  a- 1 


where  s“and  ma  are,  respectively,  the  slip  direction  and  slip  plane  normal  in  the  current 
configuration.  The  symmetric  part  of  Lpis  the  plastic  part  of  the  rate  of  deformation  tensor  and 
the  anti-symmetric  part  is  the  plastic  spin,  defined  respectively  by 


Nsys 

dp  =  -  ^  ya  (sa  (x)  ma  +  ma 

a-1 


N. 


sys 


=  ^Ya 


a-1 


and 


Nsys 

wp  =  -  ^  ya  (sa  ®  ma  —  ma  ®  sa )  . 

a=l 


(10) 


(U) 


2.3  Crystal  Slip  Constitutive  Relations 

The  crystal  slip  constitutive  model  relates  the  loading  on  the  slip  systems  to  the  slip  rate.  A 
traditional  approached  is  used,  and  the  model  is  similar  to  that  employed  by  Peirce,  Asaro,  and 
Needleman  (1983).  The  loading  variable,  ra,  is  constructed  to  be  work  conjugate  to  the  slip  rate 
so  that  the  plastic  dissipation  on  the  slip  systems  is  equal  to  the  plastic  dissipation  expressed  in 
terms  of  the  continuum  variables: 

NSys  NSyS 

^  Ya  Ta  ,  (12) 

a-1  a-1 

where  Pa  is  defined  in  equation  10. 

The  resolved  shear  stress,  ra,  provides  loading  on  the  slip  system,  which  results  in  a  slip  rate.  A 
simple  power  law  rate  model  is  used  here 


The  strength  of  the  slip  system,  ga,  is  typically  a  function  of  the  deformation,  and  a  simple 
power  law  form  is  assumed, 

ga  =  g%(l  + py)n .  (14) 

where  y  is  the  accumulated  slip, 


a:  dp 


ya  [  <t:  Pa] 
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In  all  of  the  simulations,  g$  =  33  MPa,  n  =  0.2,  fo  =  1  s  m  =  0.002,  and  p0  =  2.7  g/cm3. 
For  the  simulations  with  strain  hardening  in  section  3.3,  (3  =  20.  For  the  simulations  in  sections 

3.1  and  3.2,  the  strength  is  assumed  to  be  constant  by  setting  (3  =  0.  This  facilitates  focus  on  the 
new  aspects  of  the  model  without  the  difficulty  of  deconvolving  slip  continuity  effects  from 
strain  hardening  effects. 

A  complete  anisotropic  elastic  treatment  is  available  (Becker,  2004).  However,  as  with  strain 
hardening,  inclusion  of  anisotropic  elastic  effects  may  create  ambiguities  when  attributing 
behaviors  to  the  new  model  features.  Consequently,  the  elastic  constants  used  in  all  simulations 
create  an  isotropic  elastic  response.  The  shear  modulus  is  taken  to  be  27  GPa  and  the  bulk 
modulus  is  58  GPa. 

2.4  Non-local  Slip  Modeling 

The  introduction  slip  continuity  between  neighboring  finite  elements  is  the  main  focus  of  this 
work.  With  the  traditional  crystal  plasticity  model  and  the  finite  element  formulations  described 
previously,  the  deformation  of  neighboring  elements  is  only  connected  through  shared  nodal 
displacements  and  force  equilibrium  at  the  nodes.  Even  though  dislocations  associated  with  slip 
travel  through  the  material,  the  transmission  of  dislocations  from  one  element  to  another  is  not 
represented,  and  the  resulting  coordination  in  slip  deformation  is  not  captured.  The  effective 
result  is  that  all  finite  element  boundaries  are  infinite  sources  and  sinks  for  dislocations.  When 
viewed  in  terms  of  the  elastically  unloaded,  intermediate  configuration  described  in  section  2.2, 
all  of  the  discontinuity  associated  with  the  intermediate  configuration  is  located  at  the  element 
boundaries.  Gradients  of  lattice  orientation  within  elements  are  not  captured  when  constant  strain 
elements  are  used. 

The  goal  is  to  create  the  simplest  model  possible  that  captures  the  mutual  effect  of  slip  continuity 
on  neighboring  elements.  This  is  conceptually  most  straightforward  for  edge  dislocations  in  an 
idealized  2-D  geometry  where  the  directions  of  dislocation  motion  and  slip  are  coincident.  Screw 
dislocations,  in  which  the  material  motion  is  orthogonal  to  the  dislocation  motion,  do  not  exist  in 
the  2-D  crystal  representation.  These  would  create  shear  out  of  the  plane. 

2.4.1  Slip  Continuity 

The  simplest  approach  is  to  enforce  continuity  of  slip  without  changing  anything  else  in  the 
crystal  plasticity  formulation.  The  basic  picture  is  that  edge  dislocations  move  along  the  slip 
direction  from  one  element  to  another.  The  flux  of  dislocations  crossing  an  area  (one  in-plane 
dimension  and  the  out-of-plane  dimension,  w,  assumed  to  be  unity)  for  each  slip  system  is 

0  =  Pdis  vasa  '  n  >  (16) 
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where  pdis  is  the  dislocation  density  (number  of  dislocations  per  unit  in-plane  area),  va  is  the 
dislocation  velocity,  and  n  is  the  outward  normal  to  the  element  face.  Continuity  is  enforced  by 
requiring  that  the  flux  exiting  one  element  through  a  face  equal  the  flux  entering  the  neighboring 
element  through  the  same  face.  The  dislocation  density,  velocity,  and  the  Burgers  vector,  b 
(displacement  caused  by  the  passing  of  a  dislocation),  are  related  to  shearing  rate  by  Orowan’s 
equation: 

Ya  =  padis  va  b.  (17) 

For  a  shared  element  face,  and  assuming  that  the  Burgers  vector  is  constant,  the  flux  continuity 
in  the  idealized  2-D  crystal  can  be  represented  as 

(yasa)e  1  •  n  =  ( yasa)e2  ■  n  ,  (18) 

where  subscripts  el  and  e2  denote  the  two  elements  sharing  a  face  and  n  =  nel  =  —  ne2.  The 
slip  rates  and  slip  direction  are  generally  different  for  two  adjacent  elements  in  an 
inhomogeneous  deformation  field.  The  constraint  can  be  enforced  naturally  if  the  primary  flux 
variable  for  the  solution  were  associated  with  the  element  faces.  This  approach  is  pursued. 

In  an  implicit  formulation,  this  type  of  constraint  would  be  enforced  through  the  global  system  of 
equations.  Such  a  global  system  solution  is  impracticable  in  an  explicit  dynamic  approach  where 
many  time  steps  are  required  to  track  the  wave  motion. 

Here,  an  operator  split  approach  is  applied,  where  dislocation  fluxes  are  determined  on  element 
faces  during  a  first  phase,  and  the  deformation  due  to  those  fluxes  is  applied  in  the  subsequent 
phase.  This  is  facilitated  by  creating  a  two-pass  material  model.  The  first  pass  is  the  standard 
crystal  plasticity  model,  outlined  in  sections  2.2  and  2.3,  except  that  the  only  results  retained 
from  the  model  are  the  slip  rates.  The  stresses  and  updates  to  the  history  variables  are  discarded. 
The  slip  rates  from  this  first  phase  are  averaged  on  the  faces,  giving  values  denoted  as  yaf, 
where  the  superscript  f  refers  to  the  face  number  associated  with  the  element. 

The  slip  rate  for  the  element  (either  three-node  triangle  or  four-node  quadrilateral)  is  then 
determined  from  average  dislocation  flux  through  the  faces  as 

3  or  4  /3  or  4 

ya  =  ^  (yaf\sa  ■  it? |  if)  /  ^  (|sa  •  nf  \ if)  .  (19) 

where  if  is  the  length  of  the  element  face.  This  resulting  element-centered  slip  rate  is  then  used 
directly  in  a  modified  crystal  plasticity  module  to  update  the  stress  and  history  variables  for  the 
next  time  step. 

Free  surfaces,  where  dislocations  exit  without  constraint,  are  treated  as  natural  boundary 
conditions  by  this  approach.  No  special  boundary  conditions  need  to  be  applied  for  free  surfaces. 
Stiff  boundaries,  such  as  elastic  substrates  or  hard  precipitates,  and  tightly  bound  second  phases, 


7 


can  be  modeled  by  zeroing  the  slip  rate  at  those  interfaces.  Surfaces  with  intermediate  constraint 
that  may  allow  the  passage,  generation,  or  absorption  of  some  dislocations,  such  as  weak 
precipitates  or  grain  boundaries,  could  be  simulated  by  reducing  the  slip  commensurate  with  the 
resistance  of  the  boundary.  Ideally,  Robertson’s  rules  (Lee  et  al.,  1989)  for  transmission  of 
dislocations  across  grain  boundaries  could  be  implemented. 

2.4.2  Slip  Gradient  Effects 

A  finite  element  in  a  Lagrangian  crystal  plasticity  simulation  represents  a  fixed  material  volume 
with  dislocations  fluxing  into  and  out  of  the  volume  and  dislocation  generation  within  the 
volume.  With  the  slip  rates  on  faces  used  as  primary  solution  variables,  a  straightforward 
application  of  Reynolds  transport  theorem  over  the  volume  V  and  enclosing  surface  S, 

pSlsdV  =  f^~dV  +  I  pS,s  vas“  ndS  ,  (20) 

V  V  S 

can  be  used,  along  with  dislocation  nucleation  and  annihilation  rates,  to  estimate  a  lower  bound 
on  the  population  of  statistically  stored  dislocations  in  an  element  volume.  The  first  term  on  the 
right-hand  side  represents  generation  of  dislocations  within  an  element  and  the  second  term 
represents  accumulation  due  to  flux  across  the  element  boundaries.  Applying  the  divergence 
theorem,  and  recognizing  that  the  slip  direction  is  constant  throughout  the  element  volume  for 
the  chosen  element  types,  the  last  term  can  be  manipulated  as 


j  Pdis  -ndS  =  j  V  •  (p£s  vasa )  dV 
s  y 


f  Sa-V(p5is  va)  dV  = 

V 


r 

J 


b-lsa.v^a  dy  (21) 


The  third  integrand  of  equation  21  indicates  the  gradient  in  dislocation  flux  projected  along  the 
slip  direction,  and  the  last  integrand  relates  the  dislocation  flux  gradient  to  the  gradient  in  slip 
rate  through  Orowan’s  equation.  The  flux  can  be  evaluated  directly  from  the  surface  integral  as 


/ 


3  or  4 


3  or  4 


Pdis  VaSa 


n  dS 


-I 

/= i 


Pdis  VaSa 


nJ 


l'w  =  £ 


Yafi 


n 


f  lf  b 


-1 


w 


f= 1 


(22) 


where  w  is  the  out-of-plane  depth.  Taken  together,  equations  21  and  22  give  the  gradient  of 
dislocation  flux  along  a  slip  plane  in  terms  of  the  associated  slip  rates  on  the  element  faces. 


The  dislocations  stored  within  an  element  volume  are  associated  with  an  imbalance  of 
dislocations  entering  and  leaving  the  element,  which  is  related  to  the  dislocation  flux  gradients  of 
equation  2 1 .  If  occurring  on  a  single  slip  plane,  these  dislocations  create  pileups  leading  to  the 
Hall-Petch  relation.  If  the  excess  dislocations  are  on  many  parallel  slip  planes,  they  are  often 
associated  with  geometrically  necessary  dislocations  accommodating  the  crystal  lattice 
misorientation  between  subgrains.  These  subgrain  walls  serve  as  barriers  to  dislocation  motion, 
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so  both  scenarios  are  hardening  mechanisms.  The  potency  of  these  mechanisms  is  related  to  the 
number  of  dislocations  over  a  distance,  and  this  provides  a  material  length  scale.  Many  of  the 
extensions  to  the  crystal  plasticity  model  are  based  on  pile-ups  and  geometrically  necessary 
dislocations  with  the  intent  of  incorporating  a  size  scale  into  the  model. 


As  indicated  earlier,  the  stored  dislocations  result  from  the  dislocation  flux  gradients  integrated 
through  time.  Using  Orowan’s  relation,  the  dislocation  flux  gradients  are  related  to  slip  rate 
gradients,  as  shown  in  equation  21.  In  the  current  treatment,  the  content  of  stored  dislocations  is 
estimated  by  assuming  a  simple  time  integration  of  the  slip  rates  from  equations  21  and  22: 


r 


1 

V 


sa  •  Vya  dV 


3  or  4 

sa  •  Vya  dV  =  —  ^  yaf  sa 

V  L—k 

f= i 


n.f  l^w  =  Vy“. 


(23) 


The  time  integration  is  only  approximate  because  the  path  dependent  aspects  of  the  solution  are 
neglected  in  this  approach.  Vya  is  a  scalar  measure,  with  units  of  inverse  length,  of  the 
accumulated  slip  gradient  projected  along  the  slip  plane.  It  serves  as  a  surrogate  for  the 
dislocation  content,  either  pipe-ups  or  geometrically  necessary  dislocations. 


The  slip  gradient  is  multiplied  by  the  shear  modulus,  /i,  and  a  fixed  material  length  scale  factor, 
A,  as  a  modification  to  the  slip  system  strength  given  in  equation  14: 


N, 


sys 


9a  =  9q 


1  + 


ajiA  1  Vyc 


(i  +  P  y)n 


(24) 


a-1 


The  parameter  a  governs  the  strength  of  the  gradient  contribution.  For  the  simulations  where 
gradient  effects  are  included,  aA  =  0.001  pm,  and  a  =  0.0  for  the  simulations  where  the  slip 
gradient  does  not  affect  the  hardening. 


3.  Results 


The  primary  result  from  the  fiscal  year  2011  (FY 11)  effort  on  this  DRI  is  the  non-local,  explicit 
finite  element  code  running  a  modified,  operator-split,  crystal  plasticity  model,  where  dislocation 
flux  continuity  is  enforced  through  common  variables  at  element  boundaries  and  the  crystal  flow 
strength  is  modified  based  on  calculated  slip  gradients.  The  code  has  restart  capability  and  writes 
plot  files  in  XDMF  format  (XDMF,  2011)  for  viewing  with  a  standard  finite  element  post¬ 
processor.  Since  the  focus  of  the  work  is  introducing  physics  into  the  crystal  plasticity  model, 
large  parallel  runs  are  not  needed  and  the  implementation  is  serial. 


9 


3.1  Verification  of  Flux  Continuity  Implementation 

The  finite  element  framework,  the  crystal  model,  and  the  gradient  calculations  have  been  verified 
through  numerous  unit  program  tests  and  large-scale  simulations  with  known  solutions  and 
behaviors.  The  finite  element  and  crystal  model  were  checked  with  standard  evaluations,  which 
are  not  presented  here.  An  example  of  a  flux  gradient  verification  simulation  is  the  simple  shear 
deformation  shown  in  figures  2-4.  Here  the  non-hardening  single  crystal  sample  is  10  pm  tall 
and  5  pm  wide  with  an  element  size  of  0.5  pm.  The  top  surface  is  given  a  constant  velocity  of 
105  pm/ps  to  the  left  and  the  bottom  3  pm  of  the  crystal  is  held  fixed.  The  intent  is  to  provide  an 
abrupt  transition  from  no  deformation  in  the  lower  one-third  of  the  crystal  to  a  simple  shear 
deformation  in  the  upper  two-thirds.  Periodic  boundary  conditions  are  applied  across  the  5-pm 
direction,  so  the  multiple  elements  in  the  v-di  recti  on  serve  only  to  verify  that  the  periodic 
boundary  conditions  are  applied  properly. 

Figure  2  shows  slip  rates  and  the  gradient  of  the  slip  rates  for  slip  systems  1  and  3  when  the 
lattice  is  initially  rotated  15°  from  the  reference  orientation.  Slip  continuity  is  not  enforced  in  this 
simulation,  and  it  provides  a  baseline  result  using  the  traditional  crystal  plasticity  model. 
Focusing  on  slip  system  1,  the  slip  rates  are  0.007829  ps  1  in  the  sheared  region  and  zero  in  the 
non-sheared  region  below  it  (figure  2a).  The  slip  rate  at  the  interface  between  the  two  regions  is 
half  of  this  value.  Multiplying  this  by  the  slip  direction  contribution  across  the  interface  (sin  15), 
and  dividing  by  the  grid  spacing  (0.5  pm),  gives  a  slip  rate  gradient  of  0.002026  ps  '-pm  '.  This 
agrees  with  computed  gradient  in  the  elements  below  the  interface  to  all  four  digits  (figure  2b). 
The  crystal  lattice  has  distorted  and  rotated  in  the  sheared  elements  above  the  interface  to  an 
orientation  of  15.88°.  The  analytically  computed  slip  gradient  for  this  orientation  is 
0.002142  ps  '-pm  '.  This  is  within  0.4  %  of  the  numerically  calculated  value,  which  also 
includes  effects  of  crystal  lattice  shear  on  the  slip  direction  that  are  not  included  in  the 
analytically  computed  gradient.  A  similar  analysis  can  be  performed  for  slip  system  3,  where 
similar  accuracy  is  attained. 

Figure  3  shows  the  result  from  a  similar  calculation  but  with  the  slip  continuity  enforced  across 
element  boundaries.  The  light  blue  element  in  figure  3a  and  the  one  immediately  below  it 
indicate  a  monotonic  slip  gradient  for  slip  system  1 .  The  gradient  is  quantified  in  figure  3b, 
where  the  colors  indicate  that  the  slip  changes  monotonically  away  from  the  deformation 
discontinuity.  Slip  system  3,  on  the  other  hand,  shows  a  non-monotonic  slip  rate  pattern  near  the 
interface  (figure  3c).  This  is  more  clearly  seen  in  figure  3d  as  the  color  changes  jump  to  extreme 
values  on  either  side  of  the  elements  at  a  y-coordinate  of  2.75  pm.  This  lack  of  monotonicity  is  a 
subject  for  later  discussion. 
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Figure  2.  Simple  shear  of  a  crystal  rotated  15°  with  the  bottom  3  jam  held  fixed.  Slip  continuity 

is  not  enforced,  (a)  Slip  rate  on  system  1,  (b)  slip  gradient  on  system  1,  (c)  slip  rate  on  system  3, 
and  (d)  slip  gradient  on  system  3. 
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Figure  3.  Simple  shear  of  a  crystal  rotated  15°  with  the  bottom  3  mm  held  fixed.  Slip  continuity 

is  enforced,  (a)  Slip  rate  on  system  1,  (b)  slip  gradient  on  system  1,  (c)  slip  rate  on  system  3,  and 
(d)  slip  gradient  on  system  3. 
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Figure  4  shows  results  from  a  similar  simple  shear  calculation  with  the  crystal  lattice  in  the 
reference  orientation  such  that  one  of  the  slip  directions  is  horizontal,  parallel  to  the  shearing 
direction.  Here  there  is  only  slip  on  system  1,  so  only  results  from  slip  system  1  are  presented. 
Due  to  the  orientation  of  the  crystal  with  respect  to  the  grid,  dislocations  would  not  cross  the 
horizontal  element  boundaries,  and  the  slip  continuity  algorithm  should  not  couple  the 
deformation  of  neighboring  elements.  Figure  4a  shows  no  indication  of  slip  transmission,  and 
this  is  confirmed  by  the  gradient  plot  in  figure  4b.  This  verifies  that  the  simulations  are  only 
enforcing  dislocation  flux  continuity  when  dislocations  move  across  the  element  interfaces. 


Figure  4.  Simple  shear  of  a  crystal  rotated  0°  with  the  bottom  3  mm  held  fixed.  Slip  continuity  is 
enforced,  (a)  Slip  rate  on  system  1  and  (b)  slip  gradient  on  system  1 . 

3.2  Example  of  Flux  Continuity  Effects  in  a  Complex  Deformation  Field 

The  behavior  of  the  numerical  integration,  particularly  the  operator  split  on  the  crystal  model, 
was  evaluated  using  a  punch  simulation  where  the  gradients  are  severe.  A  single  crystal,  40  pm 
wide  and  20  pm  high,  was  deformed  as  if  indented  by  a  flat  punch  over  the  center  third  of  the 
crystal.  The  punch  displacement  rate  is  105  pm/ps.  The  grid  spacing  for  the  quadrilateral  mesh  is 
0.25  pm  in  both  x-  and  y-directions.  The  bottom  surface  is  prevented  from  vertical  motion  and 
the  lateral  surfaces  are  traction  free.  The  initial  orientation  of  the  crystal  lattice  is  15°  from  the 
reference  orientation  so  that  the  deformation  field  is  asymmetric. 

Plots  of  slip  rate,  on  a  logarithmic  scale,  are  shown  in  figure  5  for  the  standard  crystal  plasticity 
model  and  the  model  with  dislocation  flux  continuity  enforced.  As  expected,  enforcing 
dislocation  flux  continuity  smoothes  the  fields  along  the  highly  strained  bands  and  below  the 
punch  (figure  5b),  but  it  does  not  appear  to  be  excessively  dispersive.  It  also  does  not  appear  to 
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diminish  the  peak  slip  rates  significantly  when  comparing  the  sheared  bands  emanating  from 
either  side  of  the  punch. 
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Figure  5.  Slip  rate  contours  for  single  crystal  initially  orientated  15°  (a)  without  slip  continuity 
enforced  and  (b)  with  slip  continuity  enforced. 

Contours  of  crystal  lattice  rotation,  with  and  without  dislocation  flux  continuity  enforced,  are 
shown  in  figure  6.  The  enforcement  does  little  to  the  positive  lattice  rotations  on  the  left  side  of 
the  punch  where  the  rotation  is  focused  in  a  narrow  band.  However,  the  peak  negative  rotation  is 
significantly  diminished  by  enforcing  slip  continuity,  and  the  rotation  field  below  the  punch  is 
also  more  diffuse.  This  is  consistent  with  the  spreading  of  the  slip  rate  contours  in  figure  5  and  is 
expected  from  the  additional  constraints  imposed  on  the  crystal  deformation. 
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Figure  6.  Crystal  lattice  rotation  contours  for  single  crystal  initially  orientated  15°  (a)  without  slip 
continuity  enforced  and  (b)  with  slip  continuity  enforced. 

The  punch  load  for  the  simulation  with  flux  continuity  imposed  is  approximately  0.5%  lower 
than  the  punch  load  from  the  standard  crystal  plastic  simulation.  Since  the  slip  system  strength  is 
prescribed  to  be  constant  in  these  calculations,  the  lower  load  implies  that  the  crystal  lattice  is  in 
a  slightly  softer  orientation  when  dislocation  flux  continuity  is  enforced.  Understanding  this  will 
require  further  exploration. 

3.3  Exploration  of  Gradient  Effects  on  Length  Scale 

A  simple  shear  simulation  on  a  thin  film  is  run  as  a  preliminary  assessment  of  the  effects  of 
increasing  slip  resistance  due  to  gradients  in  dislocation  density.  Strain  hardening,  J3=  0.5,  is 
included  to  smooth  inhomogeneities  in  the  deformation  field  triggered  by  noise  in  the  explicit 
solution.  Three  sets  of  simulations  are  run:  the  baseline  crystal  plasticity  model,  enforcing 
dislocation  flux  continuity,  and  inclusion  of  gradient  hardening  in  addition  to  enforcing 
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dislocation  flux  continuity.  The  lattice  orientation  is  initially  15°,  and  periodic  boundary 
conditions  are  applied  so  that  only  a  single  column  of  elements  is  needed.  The  film  is 
sandwiched  between,  and  perfectly  bonded  to,  two  stiff,  horizontal  dies  that  prohibit  dislocation 
transmission.  This  creates  a  slip  gradient  that  changes  as  the  thickness  of  the  film  is  varied. 
Simulations  are  run  for  a  10-pm-thick  and  a  1  -prn-thick  film  to  ascertain  the  size  scale  effect. 

The  slip  rates  are  plotted  in  figure  7  for  the  10-pm-thick  film  and  in  figure  8  for  the  1  -pun-thick 
film  at  a  shear  strain  of  0.05.  As  expected,  the  standard  crystal  plasticity  model  (figures  7a  and 
8a)  is  strained  uniformly  through  the  thickness  and  the  behavior  is  the  same  for  the  two  film 
thicknesses.  Since  the  model  shears  uniformly,  the  sides  of  the  model  are  linear.  With  dislocation 
flux  continuity  enforced  at  element  boundaries  and  zero  flux  enforced  at  the  film  boundaries 
(figures  7b  and  8b),  the  slip  is  not  uniform  through  the  thickness.  The  slip  rate  in  the  elements 
near  the  surface  is  reduced  due  to  the  flux  boundary  condition,  and  the  slip  rate  is  highest  in  the 
second  element  from  the  surfaces.  While  there  is  some  difference  in  slip  rate  with  film  thickness, 
it  is  not  significant  compared  to  the  strength  of  the  gradient  near  the  surfaces.  The  sides  of  the 
model  are  not  straight.  There  is  a  kink  near  the  surfaces,  where  the  slip  rate  is  highest,  and  the 
center  portion  appears  linear. 

The  addition  of  hardening  associated  with  slip  gradients  modifies  the  deformation  patterns  in 
figures  7c  and  8c.  The  total  slip  rate  varies  continuously  through  the  thickness,  and  the  specimen 
profile  is  sigmoidal.  The  slip  rate  is  relatively  uniform  in  the  center  of  the  1 0- pm -thick  film,  and 
the  gradient  is  concentrated  at  the  surfaces.  The  gradient  is  spread  over  the  entire  thickness  of  the 
1- pm- thick  film.  This  is  also  reflected  in  the  greater  curvature  of  the  thinner  film. 
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Figure  7.  Slip  rate  sum  for  10-jam-thick  film  simulations:  (a)  traditional  crystal  plasticity, 
(b)  with  slip  continuity  enforced,  and  (c)  with  slip  continuity  and  strengthening 
from  slip  gradients. 
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Figure  8.  Slip  rate  sum  for  1-pm-thick  film  simulations:  (a)  traditional  crystal  plasticity,  (b)  with  slip 
continuity  enforced,  and  (c)  with  slip  continuity  and  strengthening  from  slip  gradients. 
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Slip  gradients,  at  a  shear  strain  of  0.05,  are  shown  in  figure  9  for  the  10-pm-thick  film  and  in 
figure  10  for  the  1  -pm -thick  film.  Gradients  are  presented  for  slip  systems  1  and  3  for  the 
calculations  with  dislocation  flux  continuity  enforced  and  the  calculations  with  flux  enforced 
plus  gradient  effects  on  hardening.  The  gradients  from  the  standard  crystal  model  are  identically 
zero  and  are  not  shown.  The  gradients  reflect  the  slip  distributions  given  in  figures  7  and  8.  Just 
enforcing  dislocation  flux  continuity  affects  only  the  elements  near  the  surface,  and  there  is  no 
appreciable  size  scale  effect.  When  the  slip  gradient  modifies  the  strain  hardening,  the  gradient 
on  slip  system  3  is  significantly  different  for  the  two  film  thicknesses.  The  gradient  is  smooth 
and  monotonic  for  the  thinner  film  (figure  10c)  but  not  for  the  thicker  film  (figure  9c). 


Figure  9.  Slip  gradient  on  systems  1  and  3  for  10-pm-thick  film  simulations:  (a)  and  (b)  with  slip  continuity 
enforced  and  (c)  and  (d)  with  slip  continuity  and  strengthening  from  slip  gradients. 
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Figure  10.  Slip  gradient  on  systems  1  and  3  for  1-pm-thick  film  simulations:  (a)  and  (b)  with  slip  continuity 
enforced  and  (c)  and  (d)  with  slip  continuity  and  strengthening  from  slip  gradients. 

The  effect  of  dislocation  flux  continuity  and  gradient  hardening  is  illustrated  by  the  stress-strain 
curves  in  figure  1 1 .  The  response  of  the  standard  crystal  model  is  the  same  regardless  of  film 
thickness.  The  enforcement  of  dislocation  flux  continuity  between  the  elements  modifies  the 
response  slightly,  but  not  enough  to  create  a  distinct  line  on  this  plot.  Thus,  curves  from  four  of 
the  six  calculations  appear  coincident  in  figure  1 1 .  The  effects  of  including  gradient  hardening 
are  evident.  The  response  of  the  10-pm-thick  film  departs  perceptibly  from  the  baseline  model, 
and  the  stress  in  the  thinner,  1-pm-thick,  film  is  substantially  greater.  It  should  be  noted  that  if 
the  crystal  is  initially  in  the  reference  orientation,  with  the  crystal  slip  parallel  to  the  boundaries, 
the  dislocation  flux  across  the  horizontal  element  faces  is  zero.  Consequently,  there  would  be  no 
flux  gradient  across  the  elements  and  no  size  effect. 
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Figure  11.  Stress  strain  curves  from  1 -jim-thick  and  10-jam-thick  films  with 
standard  crystal  plasticity,  with  dislocation  flux  enforced  and 
with  flux  enforced  plus  gradient  strengthening. 


4.  Discussion  and  Future  Directions 


The  crossed-triangle  elements  were  implemented  as  a  way  to  align  element  boundaries  with  slip 
planes  in  an  effort  to  capture  sharper  gradients  in  the  deformation  field.  Unfortunately,  the 
crossed  triangles  still  appear  to  give  stiffer  response  than  the  constant  stress  quadrilaterals  in 
most  circumstances  (results  not  shown).  For  large  meshes  using  quadrilateral  elements,  the  ratio 
of  number  of  degrees  of  freedom  to  the  number  of  incompressibility  constraints  enforced  by  the 
elements  is  2.  The  crossed  triangles,  through  a  geometric  serendipity  (Nagtegaal  et  al.,  1974),  are 
not  over-constrained  like  unstructured  triangular  meshes.  However,  the  ratio  of  degrees  of 
freedom  to  incompressibility  constraints  is  only  4/3.  Consequently,  the  overall  response  is 
somewhat  stiffer  than  the  quadrilaterals  even  though  the  slip  aligns  with  element  boundaries. 
These  elements  will  continue  to  be  evaluated  as  advancements  are  made  to  the  code,  since 
having  element  boundaries  aligned  with  the  slip  planes  could  still  prove  useful. 

Enforcing  dislocation  flux  continuity  by  the  operator  split  approach  described  in  section  2.4 
appears  appropriate  for  smoothly  varying  slip  fields,  as  in  the  punch  problem,  but  the  abrupt 
discontinuity  evidenced  in  the  simple  shear  examples  is  problematic.  The  oscillatory  behavior 
has  the  appearance  of  an  hourglass  mode.  But  it  could  also  indicate  that  an  upwinding  scheme  is 
needed  to  deal  with  the  sharp  advecting  gradients.  The  nature  of  the  oscillatory  behavior  needs  to 
be  explored  further  before  proposing  a  remedy.  A  susceptibility  to  hourglassing  could  be 
inherent  in  the  chosen  variable  placement  on  the  element  combined  with  a  lack  of  constraint  on 
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the  slip  gradients  at  element  faces.  Physical  stabilization  of  the  hourglass  modes  would  involve 
penalizing  slip  gradients  from  element  to  element  by  introducing  additional  constraints  at 
element  faces.  Upwinding  would  involve  consideration  of  longer-range  gradients  or,  perhaps, 
borrowing  ideas  from  the  Streamline-Upwind  Petrov-Galerkin  methods  used  in  finite  element 
simulations  of  fluids.  In  either  case,  further  development  of  the  numerical  approach  will  be 
needed.  This  could  prove  challenging. 

Enforcing  slip  continuity  does  smooth  the  deformation  field,  as  anticipated.  However,  slip 
continuity  in  the  crystal  plasticity  model  does  not,  by  itself,  appear  to  be  sufficient  to  induce  a 
proper  length  scale  dependence  in  the  simulations.  Thus,  the  simple  enhancement,  while 
physically  appropriate,  is  likely  not  adequate.  Extending  this  framework  to  include  gradient 
terms  in  the  hardening  relation  to  represent  dislocation  pileups  and  geometrically  necessary 
dislocations  does  provide  a  length  scale,  consistent  with  results  from  numerous  prior  studies 
(e.g.,  Arsenlis  et  al.,  2004;  Gerken  and  Dawson,  2008;  Gurtin  et  al.,  2007;  Mayeur  et  al.,  2011). 
However,  even  with  gradient  hardening,  the  initial  yield  strength  of  the  single  crystal,  indicated 
in  figure  1 1 ,  is  unaffected.  This  suggests  that  some  important  physical  mechanism  is  still 
missing. 

The  current  approach  is  distinct  from  the  past  work  in  the  primary  variables  used  to  calculate  the 
gradients  and  the  implications  for  boundary  conditions.  The  proposed  method  uses  fundamental 
quantities  already  existing  in  the  solution,  i.e.,  crystal  slip  rates,  and  physically  intuitive 
boundary  conditions  can  be  applied.  Other  approaches  introduce  new  degrees  of  freedom,  such 
as  dislocation  density  or  geometrically  necessary  dislocations,  which  have  no  obvious  boundary 
prescription.  This  has  been  an  impediment  to  their  use  that  could  be  ameliorated  by  the  current 
approach. 
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6.  Transitions 


It  is  envisioned  that  a  continuum  crystal  plasticity  model  that  enforces  dislocation  flux  continuity 
between  adjacent  elements  will  introduce  a  length  scale  into  grain-scale  continuum  simulations, 
thereby  enhancing  multiscale  modeling  capabilities.  The  model  could  be  used  in  grain-scale 
deformation  simulation  of  metals,  ceramics,  and  explosive  materials  to  represent  constraints  of 
grain  boundaries  and  surfaces  more  accurately.  The  capability  aligns  with  crystal  plasticity 
efforts  envisioned  for  the  Materials  in  Extreme  Dynamics  Environment  Cooperative  Research 
Alliance,  anticipated  to  begin  in  fiscal  year  2012  (FY12),  which  will  serve  as  a  venue  for 
continued  research  both  across  the  Weapons  and  Materials  Research  Directorate  (WMRD)  and 
with  the  Alliance. 
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